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THE  CHEMICAL  EQUILIBRIUM  PROBLEM 


I.  Formulation 

Consider  a  beaker,  Into  which  we  will  put  known  amounts  of 

various  atoms,  ions,  radicals,  molecules,  or  other  chemical  entities. 

These  entities,  called  'inputs'  will  have  the  property  that  no  combination 

of  more  than  one  can  react  together  to  form  other  inputs.  Thus  if  we  were 

to  load  the  beaker  with  H+  ion  and  OH  ion,  we  could  not  also  use  ^0 

as  an  input,  i.e.  We  assume  that  the  entities  when  expressed  as  vectors 

are  linearly  independent.  We  order  the  inputs  i  *  l,2,...,m  and  enter 

th  th 

the  amount  in  moles  of  the  i  input  into  the  i  component  of  a  vector  b. 

These  inputs  may  react  with  one  another  in  fixed  proportions 

to  form  various  chemical  species.  If  we  number  the  species  j  ■  1,2 . n 

then  species  j  may  be  represented  by  a  vector  with  m  compenents, 

th  th 

specifying  in  its  1  component  how  moles  of  the  i  input  are  consumed 
in  the  reaction  which  forms  one  mole  of  species  j  .  Then,  if  x^  *  the 
number  of  moles  of  species  j  in  the  solution,  conservation  of  mass  demands 
that  the  following  vector  equation  be  satisfied: 

I  P,x,  -  b  (1.1) 

j-1  J  J 

If  we  define: 

n 

x  -  ][  x  (1.2) 

j-1  J 

■^/x  is  the  concentration,  in  mole  fractions,  of  the  species. 


I 


then 


There  is  a  function,  called  the  Gibbs  Free-Energy  function,  which 


expresses  the  total  electro-chemical  potential  of  the  solution.  This 
function  in  the  single  compartmented  case  is  proportional  to  z  ;  where 


n 

“  I  xj  cj  + 

j-1  J  J 


X1  - 
log  J/x 


(1.3)' 


1.  Theory  tells  us  that  for  an  ideal  chemical  solution,  the  partial 
molar  electro-chemical  potential  of  a  species  j  takes  the  form: 

Uj  =  Mj  +  RT  log  Cj  +  Zj  FX 

where  Cj  is  the  molar  concentration  of  species  j  ,  Zj  the  electric 
charge  on  each  molecule  of  species  j  ,  and  X  the  electrical  potential. 
R,  T  and  F  have  their  usual  meanings. 

On  the  range  of  concentrations  for  which  ideality  is  a  good 
approximation  to  the  real  world,  there  is  a  constant  a  ,  approximately 
the  same  for  each  species,  such  that:  Cj  =  a  ^/x  (to  be  read 
'approximately  equal  to  ' ) .  Clearly  then, 

/**?  +  *!«  \  ^  _ 

Mj  -  RT  1 -J - ^ -  +  log  a  j  +  log  J/x 


We  let: 


V,  +  z .  FX 


+  log  a 


Then  the  actual  Gibbs  free-energy  function  becomes 


l  RT  x j  Cj  +  log  Xj/x 


Thus  z  is  indeed  proportional  to  the  Gibbs  free-energy. 

Also,  our  constants  Cj  include  any  electrical  potential 
impressed  from  outside  upon  the  compartment.  The  number  is  the 

standard  partial  molar  free-energy,  and  can  be  found  in  the  literature. 
The  constant  a  converts  mole  fractions  to  molar  concentrations.  In 
dilute  aqueous  solutions  a  will  be  approximately  55.5  moles/liter. 


Willard  Gibbs  showed  that  our  beaker  will  be  at  equilibrium  when  the 

chemical  solution  achieves  the  composition  x  -  (x^ . x^)  and  Min  z 

satisfying: 


n 

l  xj 
j-1  J 

X1  - 

c.  +  log  / X 

I  J 

■ 

z (x) (Min) 

n 

I  Pj 

j-1  J 

XJ 

“ 

b  (1.4) 

XJ 

> 

0  *  j  “  1*2, . • • i 

n 

where 

X 

“  l 

j-1  J 

Methods  for  finding  equilibrium  solutions  to  chemical 
problems  based  on  this  formulation  have  had  much  success.  R.  J.  Clasen  [1] 
at  Rand  Corporation,  for  example,  has  devised  several  procedures  for  solving 
(1.4).  However,  all  such  methods  have  so  far  run  into  trouble  whenever 
degeneracy  has  occurred  during  the  computational  procedure,  i.e.  whenever 
some  Xj  becomes  zero  or  nearly  zero.  This  paper  will  propose  a  method  > 
not  at  all  discommoded  by  degeneracy. 


II.  Method  of  solution 

Suppose  we  had  an  initial  feasible  solution  x.  ■  «... »xn) 

to  (1.4)  with  the  property  that  each  component  x^  of  x  were  strictly 
positive.*  Then  by  following  the  stops  given  below,  one  can  successively 
improve  the  solution,  and  ultimately  find  a  solution  as  close  to  the 
optimal  solution  as  desired. 


*  Before  the  method  in  this  Section  can  be  applied  it  is  necessary  to  have 
at  hand  an  initial  solution  with  all  x,  >  0  .  How  to  find  such  a  solution 
when  it  exists  is  discussed  in  AppendixJl.  If  for  certain  j,  x  ■  0  must 
hold,  Appendix  1  gives  an  algorithm  for  determining  which  x.  must  be  dropped 
so  that  the  above  algorithm  can  be  initiated  on  a  subset  of  J  the  variables. 
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Step  1: 


Letting  9j  -  £  ,  j  -  1,2, ... ,n. 


and 


0  -  =  , 
x  * 


find  the  unique  y  satisfying  the  quadratic  program: 
Min  z: 


Find  y^  >_  0  , 


2  l  ej  +  I  -  1  -  log  (ye)l  y,  - 


1  Vi  ■  b 


This  is  equivalent  to  solving 


C1  -  1 

c2-J 

c  -  1 
n 


-b 


\  “ 
A/_ 
!°g  B  'e 

V 

log  'e 

n/_ 
iog  'e 


"  * m 

°1 

- 

?2 

on 

0 

, 

(2.1) 


yj  -  0  *  1  0  . 

and  n 


l  yj 

J-l  J  J 


j  "  l,2,i. . ,n. 


Step  2: 


Form  the  weighted  average  solution  y  , 


u  " 


Ax  +  (1  -  A)  y 


(2.2) 


where  X  is  chosen  to  minimize  the  value  of  z(u(X))  for  all  values 
of  X  for  which  (2.2)  is  non-negative.  Notice  from  (2.1)  that  y  is 
a  feasible  solution  to  (1.4),  and  since  x  is  feasible,  so  must  y  be 
feasible. 


X1  "  y1 

Now  compute  °*  1  +  J  x  J 


J 


If  each  | n j |  is  small 


enough  then  y  is  a  close  approximation  to  the  optimal  solution.  If 

not,  go  to  Step  3.  We  can  decide  if  each  |n  |  is  small  enough  by  the  following. 


It  will  be  shown  that  x  ■  (x^,...,xn)  is  an  optimal  solution 
to  problem  (1.4), 


Min  z(x)  ■  £  x 


j 


Cj  +  log  Vx 


s.t. 


1  Yj 


=  b 

>  0 


if  and  only  if  there  is  a  vector  II  satisfying: 

o 

.“O  o 

0  j  =  1 9  •  •  •  »n  • 


c j  +  log  ^/x°  -  n°Pj 


We  have  assumed  here  that  x^  >  0  ,  j  ■  l,2,...,n. 

Let  x^  *  (xj,...,x^)  be  the  optimal  solution  to  the  similar  problem. 


Min  w(x)  *  £  x 


j 


+  log 


s.t. 


*  Yj 


'j 


«=  b 

>  0 


Clearly,  the  feasible  solution  x  of  both  problems  are  the  same.  Further¬ 
more,  if  x  is  feasible,  then: 
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Min  w(x)  =  I  xj  |  c j  +  1°8  J 

s.t.  I  PjXj  -  b 

*J  -  0 

Hence,  by  the  foregoing: 

|z(x)  -  z(x°)  I  <.  B  , 

where  B  is  any  bound  which  satisfies: 

2  I  'j kj  -  cj I  -  2  l  tjlbjl  i' 

for  all  feasible  solutions  t  =  (t, . t  )  . 

1  n 

Such  bounds  are  not  usually  hard  to  find.  For  example,  one 
can  usually  determine  a  bound  t^  on  t  =  Jt  ;  by  inspection.  If  the 
inputs  may  combine  together  but  in  no  species  does  an  input  split  into  two 
or  more  parts,  then  the  sum  of  the  inputs  is  such  a  bound  tm  on  t  . 

Then  if  n  ”  max  I  n.  I  ,  we  can  let: 

m  }  i 

B  *  2tmr*m 


Alternately,  a  separate  upper  bound  on  each  species  may  be 


found  by: 


mj  ■  “jln  {  1/ai3  l  biaij  - 0  •  aij  +  °} 


,  th 


00  if  no  such  b^  and  a^ 


where  a^  is  the  i  component  of  the  vector  .  If  all  these  numbers 

m.  are  finite,  then  we  can  set 
J 

B  “  2  l  mj  |  Hj  | 

Whatever  B  is  finally  decided  upon,  the  following  will  be  true: 

z(x°)  >  z(x)  -  B 


Hence  we  will  be  able  to  estimate  how  close  the  solution  u  -  Ax+(1-A)y 
is  to  the  optimum. 


Step  3. 


In  theory,  since  x  is  strictly  positive,  then  u  will  also 
be  strictly  positive  (see  Section  3).  If  the  computations  are  made  in 
practice,  it  may  be  that  some  component  of  u,  say  u j ,  is  so  small  as  to 
be  negligible.  In  this  case,  set  at  some  small  lower  bound  and 

continue  the  computations.  If  there  are  several  such  u ^ ,  they  can  be  made 
all  small  and  at  the  same  time  proportional  to  exp  (IIP ^  -  Cj). 

Alternatively,  column  j  could  be  deleted  from  the  problem 
entirely.  This  of  course,  means  that  the  species  will  not  appear  at 

all  in  the  final  solution,  and  so  it  is  treated  as  if  it  were  zero  in  the  final 
solution  to  the  original  problem. 

Step  4. 

Let  the  strictly  positive  solution  u,  modified  as  in  step  3, 
take  the  place  of  x  in  step  1  and  in  (2.1). 

One  will  repeatedly  cycle  through  these  four  steps  until  the 
convergence  criterion  of  step  2  is  satisfied.  At  that  point,  the  last 
solution  found  will  closely  approximate  the  actual  optimal  solution  to  (1.4). 

Ill,  Derivation: 

The  chemical  equilibrium  problem  (1.4  or  3.1  below)  is 
difficult  to  solve  because  its  objective,  the  function  z(x),  is  non-linear. 
Linear  problems  being  easy  to  solve,  we  will  attempt  to  replace  this  problem 
with  a  linear  one.  We  do  this  in  two  stages:  first,  we  find  an  approximating 
problem  with  a  non-linear  separable  objective  which  in  turn  is  approximated  by 
a  quadratic  program  with  a  convex  separable  objective.  The  first  step  is 
accomplished  by  theorem  1. 
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Theorem  1 

Consider  the  following  two 

problems : 

Min  z 

n 

=  ^  Xj 
j-1  J 

cj  + 

log 

X1  -  1 

/  x 

n 

s.t. 

l  Pjxj 
j-1  2  2 

=  b 

(3.1) 

Xj 

1  0 

n 

where 

*  =  l  xj 
j-1  2 

and: 

Min  w 

* 

V 

K  + 

log  ) 

s.t. 

I  Vj 

=  b 

(3.2) 

x . 

J 

_>  0 

Then  there  is  a  particular  value  of  K  for  which  the  optimal  solution 

x°  =  (x°,...,x°)  to  (3.2)  is  also  optimal  for  (3.1). 
l  n 


Proof :  Look  at  the  optimality  conditions  for  the  two  problems:  x 

is  optimal  for  (3.2)  if  there  exists  a  vector  11°  =  (n°,...,n°)  suc^  that 
(x°,  n°)  satisfies: 


a  w(x°)  -  n°  l  p.  x.  -  b 

, _ \k  I]  ■=  Cj  +  1  +  K  + 


log  Xj  -  n  Pj 


2.  0  j  =  1 . n 

(3.3) 

=  0  if  xj  >  0 

9 


Similarly,  it  is  optimal  for  (3.1)  if  there  is  a  -  (n  such 

that  (x°,  n1)  satisfies: 


( 


I  2.  0  j  -  1,. . .  ,n  (3.4) 

-  0  if  x°  >  0 

Now  let : 

K  =  -1  -  log  x°  (3.5) 

Substitution  of  (3.5)  into  (3.3)  shows  that  (x°,  11°),  for  this  value  of 
K  ,  satisfies  the  optimality  conditions  (3.4)  of  problem  (3.1).  Of  course 
this  also  shows  that: 

n1  -  n°  (3.6) 

QED. 

If  we  did  happen  to  know  the  value  of  x°  ,  then  theorem  1 
states  that  the  solution  to  problem  (3.2),  with  K  as  in  (3.5)  would  be  the 
desired  solution  to  (3.1).  In  addition,  (3.2)  is  separable,  as  we  wished. 

Unfortunately,  x°  is  not  known.  Instead,  we  will  use  the 
value  of  x  in  the  current  solution  to  (3.1)  to  approximate  K  in  (3.5)  and 
(3.2). 

If  our  current  feasible  (but  non-optimal)  solution  to  (3.1) 
is  x  =  (X^,...,#^)  ,  so  that  x  E  \  x^  ,  we  let: 
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(3.7) 


dj  «  Cj  -  1  -  log  x 
and  then  consider  the  problem: 

Min  w  =  l  Xj  (d ^  +  log  x^ ) 
s.t.  J  PjXj  *  b  (3.8) 

XJ  i  0 

The  approximating  quadratic  program  could  be  obtained  by 
expanding  u  log  u  into  a  Taylor  series  about  u  ■  .  However,  it  is 

also  possible  to  "linearize"  (3.8)  by  replacing  it  with  an  equivalent 
generalized  linear  program: 


Min  <p  (x,  0)  ■ 

1  Wj 

“  log 

V 

s.t. 

I  Vj 

K 

b  :  n 

Vj 

< 

1  :  <-yj> 

xj 

> 

0 

where 


ej  -  0 


is  a  variable  that  can  be  chosen  independently  of 


It  is  easy  to  show  that  from  any  feasible  solution  of  either 
(3.8)  or  (3.9),  a  feasible  solution  to  the  other  can  be  found.  Furthermore, 
the  same  holds  true  for  optimal  solutions. 


A  moments  thought  will  convince  one  that  if  (x,0)  is  a 
solution  to  (3.9)  which  does  not  satisfy: 


Vj 


=  1 


(3.10) 
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The  II  and  -y^  appearing  to  the  right  of  the  constraints 
of  (3.9)  are  Lagrange  Multipliers.  The  multipliers  H  corresponding  to 
equality  constraints,  are  unrestricted  in  sign,  but  the  -y  corresponding 
to  inequalities  of  the  sort  <_  must  satisfy: 


y3  -  0 


We  will  assume  that  the  current  solution  (x  ,  0 ) 


(3.11) 


I 


for  which  the  inequality  is  strict  for  at  least  one  j  . 


If  we  are  trying  to  satisfy  (3.14),  we  look  for  that  value 
of  0  which  minimizes  f  (0 j )  .  Since  f(0j)  strictly  convex,  it 
has  a  unique  minimum,  occuring  at  0^  =  0^  where 

+»i 


Thus  we  pick 


—  if  y.  >  0  ,  otherwise 
0  (some  arbitrary  large  number) 


(3.15) 


So  far,  we  have  only  asserted  the  existance  of  multipliers 

A 

(  n ,  y)  satisfying  (3.11)  and  (3.12).  Conditions  (3.11)  and  (3.12)  in 
general  are  too  few  to  determine  their  values.  But  equations  (3.15)  suggest 
that  y^  be  interpreted  as  an  amount  of  species  j  ,  which  together  with 
0^  =  —  would  correspond  to  a  new  feasible  solution  to  (3.9).  Thus  we 

J  yj 

ask  that  in  addition  to  (3.11)  and  (3.12),  the  multipliers  satisfy: 


l  Vj 


(3.16) 


t 


Writing  (3.12)  and  (3.16)  in  matrix  form  gives: 


y?  —  y. 


p  p 

1  2 


d2  -  log  ex 


d,  -  log  0 


d  -  log  6 
n  6  n 


(3.17) 


If  the  matrix  [P^,  P2,...,P  ]  has  rank  m  ,  it  is  easy 
to  show  that  there  is  a  unique  solution  (Il^y)  to  (3.17),  However,  that 
solution  need  not  satisfy  y^  _>  0  .  Thus  we  relax  the  restrictions 


slightly. 


y2  —  yn 


pi  Po  -  P 

12  n 


dx  -  log  ex 


d„  -  log  0 


d  -  log  0 
n  n 


(3.18) 


V.  >  0  .  6.  >  0  .  A  -  1.2 . n 


It  can  be  shown  by  the  application  of  complementary  pivot 
theory,  or  by  noting  that  (3.18)  is  equivalent  to  a  positive  definite 
quadratic  program,  that  (3.18)  has  a  unique  solution  (H,y)  . 

A 

Theorem  2:  If  the  feasible  solution  to  (3.8)  y  4  x  found  by  solving 

(3.18)  then  y  yields  a  lower  value  of  w  than  the  current  solution 

i  m  (^. , . . .  ,51  )  . 

x  n 

Proof :  Taking  (II, y)  from  (3.18),  we  have  that: 

f<V  '  dJ  -  log  *  "VVj-  SJ10 

We  have  shown  that 
-  log  Gj  +  yj  6  j  >.  log  yj  +  1 

Now: 

°-IVi  ’  I  Vdj  -  108  6j  -  npj  +  yiV 

*  w(x)  -  nb  +  l 

since  6  *  4  .  Similarly,  summing  all  j  for  which  y.  >  0  , 

0  -  bft  -  I  yjWj  -  log  9j  -  nPj  +  yjBj) 

“  +  108  yj  "  nPj  +  ^ 

>_  w(y)  -  lib  +  £y^ 
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Thus : 


Or: 


w(x)  -  lib  +  £  y  _>  w(y)  -  nb  +  J  y^ 

a 

w(y)  <_  w(x) 


And  the  inequality  is  strict  if  there  is  any 


yj  ” 0 


such  that  6 


j 


Q.E.D. 

However,  problem  (3.8)  is  not  the  problem  we  are 
interested  in.  We  must  show: 

Theorem  3:  y  found  by  (3.18)  yields  a  lower  value  of  z  in  (3.1) 

than  does  the  current  solution  x  . 


Proof : 


Inserting  d^  from  (3.7)  into  f(0j)  » 


Evaluating 


f(0j)  *  Cj-1  -  log  x  -  log  0j  -  np  0  - 
l  Xj  f  (6j )  and  £  y ^  f (^  )  ,  we  find 

l  Xj  f  (0j )  -  z(x)  -  x  +  l  y  -  lib  -  l  Xj6j  _>  0 


I  y*  f(“  )  ■  z(y)  -  nb  +  y  log  (y/x)  <  0 
J  yj 


A  A  ^ 

!,  we  expressed  £  x.  f(0.)  and  [  y  f(-  )  in 
J  3  J  yj 


In  theorem  2, 

terms  of  w(x)  and  w(y)  .  Simple  algebra  shows  that: 


z(x)  -  z (y)  -  w (x)  -  w(y)  +  y 
=  w(*)  -  w(y)  +  yg 


X  ,  ,  X 

7  -  1  -  log  7 

y  y 


(3.19) 
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where  g(y)*y-l-  log  y  is  defined  for  y  _>  0  . 

It  is  easy  to  check  that  the  strictly  convex  functions  g(y)  has  a 
minimum  at  y  ■  1  of  g(l)  ■  0  . 


Thus: 


g(y)  1  0 


Since:  y  -  l  y^  >  0  , 

A 

and  from  theorem  2,  w(x)  -  w(y)  >  0  ,  we  have 


z(y)  <_  z (x) 

As  in  theorem  2,  the  inequality  is  strict  if  for  any  y^  >  0  ,  the 
A  X 

corresponding  0.  4  • 

J  yj 


Q  •  £  .  D  • 


So  far  then,  we  have  a  procedure  for  finding  an  improved 
solution  y  for  (3.1),  starting  from  a  non-degenerate  solution  x  . 

In  order  to  apply  the  method  again,  we  must  find  a  new  non-degenerate 
solution. 


Theorem  4:  If  y  is  a  degenerate  solution  to  (3.1),  and  x  is  a 

non-degenerate  solution  to  (3.1),  then  there  exists  X  >  0  such  that 
the  solution  y  defined  by: 

A 

yj  -  (l-X)yj  +  Xxj 

satisfies  z(y)  <  z(y)  ,  Ay  -  b  ,  n  >  0  . 

A 

(  y  is  non-degenerate  since  x  is  ,  and  \  >  0  .) 
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I 


Proof : 


Consider  z  a  function  of  A  . 


Then: 


uj  (A^, 

Cj  +  log  x*u(A) 


z(A)  -  l  u(A)  .  (A) 

j  2 

Taking  the  derivative,  and  evaluating  at  A  ■  0  (i.e.  u  ■  y)  , 


.  y  (x  _  v  n 
dA  j  Uj  V 


y*  _ 

+  log  J/y 


Since  y  is  degenerate,  we  will  suppose  y^  -  0  .  Then 

yk  - 

-  yv  >  0  ,  and  log  /y  -  -  » 


Thus : 


dz(u(X)) 

dA 


A  -  0 


dz 


Since  z  and  —  are  continuous  functions  of  A  ,  for  1  >_  A  >  0  , 


dz 


there  will  be  an  interval  (0,  6)  in  which  ~  is  negative.  Thus: 


z(u(8))  <  z (0)  ■  z (y)  . 


Q.E.D. 

Corr •  If  there  exists  a  non-degenerate  feasible  solution  to 

(3.1),  then  the  optimal  solution  will  be  non-degenerate. 

We  would  like  to  show  that  the  value  of  z  in  (3.1) 
computed  using  the  solution  u  is  strictly  lower  than  z  computed 
using  ft  . 
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Theorem  5;  Either  y  ■  x,  in  which  cane  x  is  the  optimal  solution 

A 

to  (3.1),  or  y  4  x»  in  which  case  there  is  a  solution  of  the 
form  (3.20)  which  yields  a  value  of  z  strictly  less  than  z  evaluated 
at  x  . 


Proof : 


If  y  ■  x  ,  then: 


f(6j)  "  Cj  ~  1  "  108  9j  “  log  *  "  nPj  +yj6j 


-  C  +  log  j/x  -  IIP  -  6j  , 


remembering  that  0.  ■  -r  .  Since  x  is  non-degenerate,  so  is  y  . 


Thus 


l  -  0  ->  i}  -  0  ,  J  ■  1.2 . n. 


Hence : 


- 


Cj  +  log  J/x  -  HPj  -  0  ,  j  -  1,2 . n.  (3.21) 


But  (3.21)  are  the  optimality  conditions  for  (3.1). 


Thus  if  y  ■  x  ,  then  x  is  optimal. 

If  y  4  x  i  and  y  is  degenerate,  we  have  from  theorems 
3  and  A  that  there  is  a  u  of  the  form  (3.20)  satisfying: 

A 

z(u)  <  z(y)  <,  z(x) 

A 

If  y  4  x  and  y  is  non-degenerate,  then  the  concluding 
remark  of  theorem  3  holds,  so  that 

z(y)  <  z (x) 

In  this  case  the  best  value  of  A  might  be  either  A  >  0  or  A  -  0  . 


I 


19 


Our  process,  then,  is  an  iterative  one.  We  start  with 

A  A  A 

a  non-degenerate  solution  x  -  ^ . xr)  to  (3.1)  and  find  (y,n) 

satisfying  (2.1),  which  is  repeated  here  as  (3.22); 


y2  —  yn 


\r 

c.-l  -  log  ' 0 
9  2 

c2~l  ~  log  ^6 

:  0 

n/- 

c  -1  -  log  '0 
n 


-b 


j  "  l,2,...,n 


(3.22) 


where  0,  =  — 
j  x 


j 


1 

*  * 


Problem  (3.22)  can  be  easily  derived  by  replacing  the  constants  d^  in 
(3.18)  by  their  values  given  in  equations  (3.7). 


Then  we  form  the  new  non-degenerate  solution  u  to  (3.1) 
by  equations  (3.20),  choosing  \  to  minimize  the  value  of  the  function 

A 

z(u(A)).  The  new  solution  u  replaces  x  ,  and  the  process  is  repeated. 
Each  iteration  produces  a  strict  improvement  in  z  ,  or,  if  x  is  optimal 
reproduces  the  current  solution  x  and  the  corresponding  value  of  z  . 
(See  Sec.  2) 
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The  only  remaining  question  is  whether  the  process 
converges,  and  if  it  converges  to  the  optimal  solution. 


Theorem  6  This  process  converges  to  the  optimal  solution. 

Proof :  For  any  real  chemical  problem,  the  set  of  feasible  solutions 
x  to  (3.1)  is  bounded.  Thus  the  possible  values  of  z  are  bounded.  In 
particular,  z  is  bounded  from  below,  by  z°,  the  optimum. 


Suppose  successive  applications  of  the  method  yielded  a 
sequence  x\  x^ , . . . ,xm, . . .  of  solutions  to  (3.1)  with  corresponding  values 


of  z  , 


z'1*  ,,«>  >  ...  >  ,<■»  > 


,»> 


Since  the  sequence  {zv  '}  is  monotone  decreasing  and 
bounded  from  below,  it  must  converge.  Suppose  it  converges  to  z  >  z°. 

The  sequence  {x^  }  has  elements  taken  from  a  compact  subset 

(\) 

of  n-space.  Thus  it  possesses  a  subsequence  {x  }  which  converges, 

(nk) 

and  of  course,  {z  }  also  converges,  to  z  . 

Clearly,  since  z  is  a  continuous 


<\) 

Let  lim  x  =  x 


function  of  x  ,  we  can  say  that 

l  =  z(£)  . 

And  since  i  >  z°,  k  cannot  be  optimal. 


Thus,  by  Theorem  5,  if  we  apply  our  method  to  x  ,  we  will 

A  ^ 

find  a  solution  u  to  (3.1)  such  that  z(u)  <  z(£)  . 
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Let  z(u)  =  z(x)  -  h  ,  for  some  h  >  0  . 


I  claim  however,  that  ti  is  a  continuous  function  of  St 

Let  x  be  as  in  (3.26).  Clearly  £  is  a  continuous 
function  of  St  .  And: 

u  -  A  St  +  (1-X  )  y 
o  o  J 

where  X  is  that  value  of  X  at  which: 
o 

z(u  (X))  =  Min. 

By  a  theorem  by  Dantzig.et  al  [2]  on  the  continuity  of  the  minimum  set 
of  a  function,  u  is  a  continuous  function  of  x  .  Thus  z(u)  is  a 
continuous  function  of  z(x)  . 

Thus  it  is  possible  to  pick  a  k  so  large  that  the 

nk  nk 

solution  u  derived  by  our  method  from  x  is  as  close  as  desired 

nk  \ 

to  u  ,  since  x  will  be  close  to  x  ;  and  the  value  z(u  )  will 

nk 

be  as  close  as  desired  to  z(u)  .  Let  z(u  )  <  z(fl)  +  e  ,  for  some 
small  0  <  e  <  h  .  Then: 

nk 

z(u  )  <  z(d)  +  e  <  z(x)  -  h  +  e  <  z(k)  . 
nk  nk+1  nk+1 

But  u  =  x  ,  and  z(x  )  >  z(x) ...  contradiction. 

Thus  zn  converges  to  z°  , 

Since  the  optimal  solution  x°  to  (3.1)  is  unique, 

r  n,  t  o 

ix  }  must  converge  to  x 


Q.E.D. 


22 


APPENDIX  1 


♦ 


AN  INITIAL  NON-DEGENERATE  SOLUTION 


Before  the  method  given  in  section  2  can  be  applied,  an 


initial  non-degenerate  feasible  solution  to  the  problem  must  be  found. 


That  is,  we  must  find  x  =  (x^,...,xn)  such  that: 


P.x,  =  b 

J  j 

(Al.l) 

x .  >0 

J 


To  do  this,  first  define  the  vector  Q  ,  where: 


Q  -  IPj  (A1.2) 


Then  consider  the  problem: 


Max  y 

s.t.  Qy  +  l  IjU  =  b  :  \  (A1.3) 

"j 

This  is  a  simple  linear  program,  and  therefore  easily  solved. 
Let  (y°,  u°,...,u°)  be  the  optimal  solution. 


Suppose  y  >  0 

(Al.l)  would  be: 


Then  a  strictly  positive  solution  to 


°  ,  o 

xj  =  uj  +  y 


(A1.4) 


If,  on  the  other  hand,  there  were  no  feasible  solution  to 
(A1.3),  or  if  y°  <  0  ,  then  there  would  be  no  feasible  solution  to 
(Al.l).  This  statement  is  clear  since  if  there  were  a  feasible  x  to 
(Al.l),  then  by  setting  =  x^  in  (A1.3),  we  would  have  a  feasible  solution 
with  y  =  0  to  (A1.3). 
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Finally,  suppose  Max  y  =  0  .  We  must  find  at  least  one 
pj  which  is  constrained  to  be  zero  by  the  equations  (Al.l),  (for  other¬ 
wise  we  could  reduce  all  p^  by  A  ■  Min  p^  and  replace  y°  by 
y°  +  A  ,  implying  y°  is  not  maximum.)  Now  let  A°  be  the  vector 
of  multipliers  associated  with  the  optimal  solution  (y°,  p°)  to  (A1.3). 


Pricing  out, 


1  +  A  Q  -  0 


(A1.5) 


X°Pj  1  0  *  0  and  J"1*2 . n. 


From  the  duality  theorem  A°b  =  -y°,  hence  the  equation: 


^°Q)y  +  [  (x0Pj)  j  -  -yc 


(A1.6) 


is  an  equation  which  must  be  satisfied  by  all  solutions  to  (A1.3),  Since 
y°  *  0  ,  and  A°Q  *  -1  ,  we  have: 


-y  +  l  1*^)  pj  •  o 


I  *°p, 


A  Q  -  -1,  by 


(A1.7) 


(A1.2), 


so  that  lor  at  least  one  j  we  can  say  that: 


A  P.  <  0 
J 


(A1.8) 


Letting  =  -  A  Pj  ,  we  have  from  (A1.7), 


y  +  I  =  0 


(A1.9) 


We  have  argued  that  we  should  only  consider  solutions  to  (A1.3)  which 


satisfy  y  _>  0  ,  so  that  (A1.9)  becomes: 

l  Vj  ^ 0 


(A1.10) 
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where  a ^  0  for  each  j  ,  and  >  0  for  at  least  one  j  .  That 

is,  a  positively  weighted  partial  sum  of  non-negative  variables  must 
be  less  than  or  equal  to  zero.  This  can  only  occur  when  those  variables 
are  zero. 


Thus  we  delete  those  columns  j  from  (Al.l)  for  which 
A°Pj  <  0  ,  and  again  apply  the  method.  The  process  will  continue  either 
until  we  have  found  a  reduced  problem  (A1.3)  with  y° >  0,  or  until  we  have 
deleted  all  but  m  of  the  columns  of  (Al.l).  If  the  latter  occurs,  then 


there  is  only  the  one  feasible  solution  to  the  original  problem,  and  it 


must  be  optimal. 
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APPENDIX  2 


MULTIPLE  COMPARTMENTS 


Suppose  that  our  beaker,  in  section  1,  were  divided  into 
several  different  compartments,  each  set  off  from  the  others  by  semi- 
permeable  membranes.  Some  species  would  exist  in  all  compartments, 
others  would  be  excluded  from  some  of  the  compartments. 


This  method  can  readily  be  extended  to  find  the  equilibrium 
distribution  of  the  inputs  among  the  various  species  in  the  different 
compartments.  The  mass  balance  equations,  while  still  taking  the  form: 


I  P4X- 
L  j  J 


(A2.1) 


are  by  species  into  the  various  compartments,  thus: 

\ 


l  pW  +  I  +  ...  +  l  pVJ  -  b  (A2.2) 

j-1  32  j  =  l  3  2  j-1  J  3 


where  species  j  in  compartment  i  has  vector  description  in  terms 

of  the  inputs.  Now,  however,  an  extra  equation  must  be  added  to  the 
system  (A2.2)  for  each  compartment  beyond  the  first.  This  new  equation 
will  state  that  the  net  electrical  charge  within  a  compartment  is  zero; 


l  a*  x*  =  0  ,  i  =  2,3,. ..,k 
j=l  3  3 


(A2.3) 


where  a^  is  the  charge  per  mole  of  the  species  in  compartment  i  . 
Equations  (A2.3)  can  be  appended  to  (A2.2)  by  including  the  charge  per  mole 
of  species  j  in  compartment  i  as  another  component  of  P^  . 
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An  equation  of  the  sort  (A2.3)  is  not  needed  for  one  of  the 
compartments  because  presumably  the  inputs  b  are  themselves  neutral, 
forcing  the  entire  system  to  electrical  neutrality.  If  all  but  one 
compartment  are  individually  neutral,  and  all  compartments  together  are 
neutral,  then  the  remaining  compartment  must  also  be  neutral. 


Remember  also  that  if  a  species,  say  water  (I^O) ,  is 
permitted  to  enter  more  than  one  compartment,  its  vector  representation 


q  must  be  reproduced  in  each  compartment  where  it  is  allowed.  If 
a  species  is  represented  by  its  column  P  in  only  one  compartment,  it 
will  never  occur  in  any  other. 


The  entire  multi-compartment  formulation  becomes: 
nl  /  x  1\  n2  /  x  nk  /  x 

•  -  I  c}+log  -j-  +  I  x?  Alog  ±  ♦  ...  +  l  x*  Alog  -i 

j-1  3  \  3  x  I  J-1  3  \  3  X  I  j-1  3  \  3  X 


i.t.  I  pj  x] 
L  j  i 


♦  l  * 


•••+ 


Xj  >  0  ,  i  “  l,...,nj,  i  -  1 . k.  (A2.4) 

Now,  of  course,  as  well  as  reflecting  the  relative  free  energies 
of  the  various  species  with  respect  to  the  inputs,  the  c^  must  reflect 
any  forces  applied  to  the  entire  compartment.  A  voltage  inposed  from 
outside  on  one  compartment  but  not  another,  or  a  difference  in  mechanical 
pressure  between  compartments  or  the  outside  will  alter  the  c^  . 


Once  the  problem  is  formulated,  the  method  for  solving  it 
becomes  the  same  as  in  section  2.  Now,  however,  we  let: 
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where 


x 


-i 

x 


“i 

I 

j-1 


i 


x 


J 


Then  (2.1)  becomes:  Find  (y,  n)  satisfying: 


(A2.5) 
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The  modifications  of  the  theorems  in  section  3  and  their  proofs  necessary 
to  handle  the  multi-compartment  case  are  obvious,  except  perhaps  in  the 
case  where  all  quantities  in  a  compartment  vanish. 

For  example,  consider  the  case  below: 


Compartment  1 

Compartment  2 

RHS 

H2O  Sugar 

h2o 

IO" 

1 

1 

1 

10. 

with  all  free  energy  constants  c^  -  0  .  The  second  compartment  vanishes 
in  the  optimal  solution  for  this  problem.  However,  it  is  important  to 
note  that  the  concentration  of  H^O  in  the  second  compartment  is  never 
zero.  In  any  solution  in  which  there  is  a  positive  amount  of  ^0  in 
compartment  2,  the  compartment  2  concentration  of  H2O  is  1.  Thus  in  the 
limit,  as  we  tend  toward  the  optimal  solution,  the  concentration  of  HjO 
in  compartment  2  is  still  1. 


In  general,  if  a  given  compartment  has  not  vanished,  then 
the  sum  of  the  concentrations  of  the  species  in  that  compartment  will  be 
unity.  Thus,  in  the  limit  as  the  compartment  vanishes,  the  concentration 
of  species  in  the  compartment  must  still  sum  to  unity.  The  problem  is 
only  to  apportion  this  total  concentration  of  one  among  the  species. 


Suppose,  at  some  iteration,  all  the  quantities  y^  in 
(A2.6)  either  vanish  or  become  so  small  that  they  may  be  considered  to 
have  vanished.  That  is,  all  of  compartment  1  is  gone. 
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Harking  back  to  the  optimality  conditions  to  the  original  problem,  we 
would  like  the  concentrations  of  species  in  compartment  1  to  satisfy: 


Cj  +  log  rij  -  HPJ  >  0  ,  j  =  1,2 . n^ 


(A2.7) 


where  n ^  is  the  concentration  of  the  j  1  species  in  the  first 
compartment,  and  the  n  we  use  in  the  one  we  have  available — the  one 
computed  in  (A2.6).  In  addition,  of  course,  we  demand  that: 


l  n ,  -  1 
J-l  J 


nj  -  0  »  i  1,2 . nl 


(A2.8) 


(A2.9) 


Suppose  we  let: 


where 


hj  -  exp  [npj  -  cj]/s 


s  -  l  exp[npj  -  c£] 
k-l 


(A2.10) 


(A2.ll) 


Then  clearly  (A2.8)  and  (A2.9)  are  satisfied.  Furthermore,  from  (A2.6) 


every  species  in  compartment  1  must  satisfy: 


c]  -  1  +  -  log  J 


Ve1  -  up}  > 


(A2.12) 


Since  the  compartment  has  just  vanished,  every  species  must  have  either 
decreased  in  amount  (if  it  were  above  the  minimal  allowed  level  previously) 
or  at  least  not  increased  (if  it  were  at  the  minimal  allowed  level).  Thus 
for  each  j  ,  we  have  y^e^  —  1  so  that: 


-  log  ej/e1  >_  npj  -  cj 
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Exponentiating  and  summing: 


nl  -1  nl 

1  -  l  I  I  exp[npl  -  C1]  -  s  (A2.13) 

k-1  k-1  3  3 

Thus,  substituting  (A2.10)  into  (A2.7)  and  using  (A2.13)  we  find: 

cj  +  log  i\  -  npj  -  -  log  S  _>  0 
Equations  (A2.7)  are  also  satisfied. 


The  rest  is  simple.  To  enter  a  new  iteration  of  (A2.6), 
compute  the  new  0^  by  choosing  an  appropriately  small  quantity  for 
/01  ,  and  using  as  in  (A2.10)  in  the  equation: 


01 


(A2.14) 
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